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SUMMARY 


An inviscid nonuniform axisymmetric transonic code is developed 
for applications in analysis and design. Propfan slipstream effect 
on pressure distribution for a body with and without sting is 
investigated. Results show that nonuniformity causes pressure 

f 

coefficient to be more negative and shock strength to be stronger 
and more rearward. Sting attached to a body reduces the pressure 
peak and moves the rear shock forward. Extent and Mach profile 
shapes of the nonuniformity region appear to have little effect on 
the pressure distribution. Increasing nonuniformity magnitude makes 
pressure coefficient more negative and moves the shock rearward. 

Design study is conducted with the CONMIN optimizer for an 
ellipsoid and a body with the NACA-0012 contour. For the ellipsoid, 
the general trend shows that to reduce the pressure drag, the front 
portion of the body should be thinner and the contour of the rear 
portion should be flatter than the ellipsoid. In a uniform flow of 
Mach number equal to 1.1, a reduction in pressure drag of 14 percent 
is achieved; while at a Mach number of 0.995, only 5 percent in drag 
reduction is possible. In a nonuniform flow of Mach number 0.995 to 
1.1, a drag reduction of 13 percent is obtained. For the design of 
a body with a sharp trailing edge in transonic flow with an initial 
shape given by the NACA 0012 contour, the pressure drag is reduced 
by decreasing the nose radius and increasing the thickness in the 
aft portion. Drag reduction achieved in a uniform flow of Mach 
number equal to 0.98 is 46 percent; in a nonuniform flow of Mach 
number equal to 0.98 to 0.995, 29 percent. 
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1. INTRODUCT ION 


Typical transonic axisymmetric nonuniform flows include propfan 
flow around a nacelle and a center body immersed in a jet. By 
introducing a rotation function to account for nonuniformity 
effects, a potential-like equation can be derived from the Euler 
equation, valid along a streamline. Therefore, the problem can be 
solved by revising an existing full-potential code, such as 
Reference 1. This idea was used by Brown (Ref. 2) in the transonic 
axisymmetric nozzle problem. The same formulation is presented in 
Reference 3 for an airfoil in a nonuniform flow. In both cases, a 
total velocity function is used as the primary variable. 

Optimal axisymmetric shapes have been sought experimentally by 
Whitcomb (Ref. 4) for subsonic free stream. Based on a slender body 
theory, von Karman's ogive (Ref. 5) and Sears-Haack body (Refs. 6 
and 7) can be analytically derived. Chan (Ref. 8) coupled a 
transonic small-disturbance code (Ref. 9) with a simplex optimizer 
(Ref. 10) to determine numerically optimized shapes at uniform free- 
stream Mach numbers of .98 and 1.1. However, the transonic small- 
disturbance equation is not appropriate for computation of drag for 
shapes of the blunt -nose type frequently used at transonic speeds. 
Optimal shapes in axisymmetric nonuniform transonic flow have not 
been investigated in the past. 

In this paper, a method based on disturbance potential-like 
equation is presented to solve the nonuniform, axisymmetric 
transonic problem. It is suitable for subsonic to low supersonic 
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nonuniform flow and shapes of the blunt-nose type* Optimal shapes 
with minimum pressure drag will be sought by coupling analysis with 
an optimizer (Ref. 13), using the maximum thickness and the trailing 
edge closure as constraints. Effects of different Mach number 
nonuniformity and profile shapes will also be investigated. 
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2. THEORETICAL DEVELOPMENT 


2.1 Governing Equations for Axlsyametric Nonuniform Flow 


The steady Euler equation along a streamline is (by combining 
Equations 1, 2, and 5 in Reference 3) 

2 1 -*■ I ->■ j 2 

a V • q - j q * V|q| (1) 

where a is the local speed of sound and q is the local velocity 
vector. To satisfy the surface boundary conditions exactly, the 
body-normal coordinates are used in the nose region to fit the blunt 
nose, and sheared cylindrical coordinates are used on the afterbody 
to accommodate corners such as boattails and flares (Ref. 11). For 
smooth, closed, convex bodies which are blunt on both ends, the 
transformed coordinates ( £, n) are chosen to be the usual tangential 
and normal body coordinates. In this report, the body-normal 
coordinates are used up to the first horizontal tangent; and beyond 
that point, a sheared cylindrical system is introduced. To derive 
equations in body-normal coordinates, Equation (1) is first written 
in a general curvilinear coordinate system as 


h l h 2 h 3 


'll: < h 2 h 3 u > 


9 

3x„ 


(h 3 h 1 


v) + 


9 

9x„ 


(h^w) 


1_ ,u 9_ v 9_ w_ 

2 9x^ 9x^ h^ 


g|-) |u 1 + v j + w£| 2 


( 2 ) 


where u, v, and w are the x^, and x^ components of the velocity 
vector q; and h^, and h^ are the corresponding metrics. 
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For body-normal coordinates, the metrics are 
h^ = 1 + <n 
h 2 = 1 

h 3 = r + n cos 0 


(3a) 

(3b) 

(3c) 


where 0 and k are the angle (measured counterclockwise from the axis 
of symmetry) and curvature of the reference coordinate surface; r is 
the radius from the axis; and the corresponding coordinates are 


X 1 = 5 

(4a) 

x 2 " n 

(4b) 

X 3 - ? 

(4c) 


as depicted in Figure 1. Notice that w = 0 for axisymmetric 

cases. Now Equation (2) can be expressed as follows: 

, 2 2. 1 A . N , 2 2, 

(a - u ) I u^ - uv( f + u^) + (a - v )v n 


2,k , cos0 v , 2 sin0 

+ a (— + )v + a u = 0 

V H r ' r 


(5) 


where H = h^ = 1 + kti, and subscripts denote partial 


differentiation. Define a velocity function <f> and a rotation 


function F to relate velocity components u and v as follows: 

u = i <J> + (1 + F)cos 0 (6a) 

h 5 

v - ♦ - (1 + F)sin0 (6b) 

Then Equation (5) can be reduced to a second-order partial 
differential equation in <f> with rotation function derivatives as 
forcing functions as follows: 




4 




sin0 N 1 

H h 


+ [(1 - 




COS0 




n 


+ 


sine + (1 - ^cose] -g F 
a a 


2 

- [~ cos0 + (1 - ^)sin0]F =0 (7) 

a a 

This equation is similar to the corresponding uniform flow 
disturbance potential equation with the addition of rotation 
function derivatives as forcing functions. 

In the sheared cylindrical coordinates, the velocity function <|> 
and the rotation function F ar related to velocity components u, v, 
as 


u = 1 + F + - r b '(|> n 

(8a) 

< 

II 

-e- 

(8b) 


Thus the governing equation becomes 


< x -¥♦«- 2[r b (1 -4 +i T^n 

a a a 


+ 1(1 - ~2> + ( r b) 2 (l ” ~2> + ~1T” r l>] 

a “ 


2 J 

a a 


+ lr - 'i* 1 V -4> F C 

a a 


“ - \> + - 0 
a a 


(9) 


where £ is identified with the axial coordinate, x, and n is a 
transformed radial coordinate such that ri = 0 is the body surface. 
The body shape enters Equation (9) through the first derivative 
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r^ and the second derivative r” of the local body radius, where 
primes mean differentiation with respect to x. With the two 
coordinate systems joined as described, the body surface is a 
coordinate surface where n = 0, and this simplifies the application 
of the surface flow tangency condition. It is also observed that 
Equation (9) has the same coefficients as the uniform-flow 
disturbance potential formulation again except the F^ and F^ terms. 


2.2 Equation at the Axis 


Along the axis of symmetry (the stagnation streamline) the 
limiting form of Equation (5) must be used to properly treat the 
terms involving The following symmetry conditions are used; 

<f>£ = 0 (10a) 

F^ = 0 (10b) 

and since 0 = 90° at the axis of symmetry, 


u = 0. 

The following limits are used as £ 0: 

cos 9 Ksin0 _ k 
r * Hsine = H 

u Ell 1 k 1 

? * Hsin0" = H ( H V € " % " I +5 " "2 

li n 


(10c) 


(Ha) 

(lib) 


where the subscript o denotes a quantity on the axis. Since the 
rotation function is constant along the axis, 


F n = 0 

Hence at the axis, Equation (5) becomes 


( 12 ) 
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0 


(13) 


~J + (1 ""T^o + 2 H~ ^on 

H ^ a nn o' 

o o 

Notice that the rotation function derivatives are not present in 
Equation (13) and that Equation (13) is identical to the uniform- 
flow potential formulation. 


2.3 Rotation Function 


The vorticity vector ui is defined as the curl of the velocity 
vector and can be shown to be 


a) = V x q 

"I l v 5 - 

= F + F cose)£ (14) 

h 5 n 

Its magnitude can be further linked to thermodynamic properties as 
(Ref. 2) 


sin9 

H 


F, 


+ F cos 0 = 

n 


u(l + 





p ) 
o n 


(15) 


where R = 1 is the normalized gas constant, T is the temperature, P 
is the pressure, M is the local Mach number, y is the ratio of 
specific heats, and the subscript o denotes stagnation quantities. 
Note that all variables in the above have been implicitly normalized 
with respect to the ambient velocity q a and the ambient pressure 
P a - Define the stream function as follows: 

n 

iKn) ■ i pudt (16) 

o 

where t is a dummy variable and p is the local density. 
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Since stagnation quantities are constant along a streamline, 
the stream function i|j can be used to identify the stagnation 
pressure P the stagnation density p Q , the stagnation temperature 
T 0 , and thus the stagnation speed of sound a Q . After local Mach 
numbers are calculated, the local density can be obtained by the 
isentropic relation 

1 

p = p Q /(l + 1^1 M 2 )^ (17) 

In sheared cylindrical coordinates, the vorticity vector is 


ai = (v - u )£ = -F lc 
x y' n 

So Equation (15) becomes 


F 

n 


yR 


u( 1 + 


y - i 




8I 0 


o 


ap„ 

— ) 


(18) 


(19) 


Because Equation (19) is identical to Equation (15) if 0 is 
equal to zero, the boundaries of these two coordinates are therefore 
chosen to be the first horizontal tangent on the body surface. 

Now the solution procedures to calculate the rotation function 
can be stated as follows: 

(1) For each constant £, assume the initial local density to be 
that of the undisturbed one. 

(2) Calculate the stream function using Equation (16). 

(3) Interpolate p 0 , P 0 , and T Q ; and calculate a Q , M, and p, using 
Equation (17). 

(4) Obtain T qti and P Qr) by applying cubic spline interpolation to 
T 0 (ri) and P Q (n) values. 
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(5) Repeat steps (2) - (4) until density converges. 

(6) Integrate F^ to obtain the rotation function by moving F^ term 
to the right side of Equation (15). 

The iterative process will converge in several iterations. 


2.4 Pressure Coefficients 


Along a streamline, the following form of energy equation for a 
perfect gas can be used: 


1 2 
2 



1 

2 


<&*) + 


Y - 1 


a o ( *> 

7^~T 


( 20 ) 


If the entropy is assumed to be nearly constant along a 
streamline, i.e. only weak shocks are present, the pressure 
coefficient can be derived from Equation (20) as 


C 

P 


P 


P 

00 


ref ref 


{[i + 


Y - 1 


M*U)(1 


yM 


ref 


|qj*>l 2 


i} 


( 21 ) 


2.5 Coordinate Stretching Functions 


The normal coordinate n will be stretched according to the 
following relation (Ref. 1), 

Ay 


n = 


(i - y) 


( 22 ) 
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where y is the computational coordinate which varies from zero at 

the body to one at infinity. The constant A controls the physical 

step size at the body (denoted as o) , A = n ; and for a given value 

y o 

of A, the exponent a controls the size of the last finite value of 
n* Large values of a move points farther away from the body. 

The tangential coordinate stretching to be used is a 
transformation between the physical arc length along the reference 
surface, £, and the computational coordinate, x, which varies from 
zero to one. For closed bodies the transformation is (Ref. 1) 

5- 


5 = + (x - ^>l A + B < x “ 


(23) 


where A and B are determined by specifying 5 (° denotes nose or x 

o 

= 0) and requiring that 5=5 at x - 1. These conditions give 

max 


A = (35 - 5 ) 

2 max x 

o 

B = ~ A ) 
max 


(24a) 

(24b) 


For open bodies the tangential coordinate stretching is divided into 
two regions with the physical location of the dividing point being 
x m . The stretching function for the region from the nose up to x^ 
is given by (Ref. 1) 


r . 3 5 7 

E, = a^x + a^x + a^x + a^x 


0 < x < x (25) 

m 


In the region from x ffl to infinity, the stretching function is 
(Ref. 1) 


5 - + b 

m 


(x - x )(1 - x ) 
m m 


X < X < 1 

m 


(26) 
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The coefficients in these expressions are determined by 


specifying £ , £ , and £ and requiring that £ and E, be 

m x x x xx 

o xm 

continuous at x = x^ These conditions give (Ref. 1) 


a l ? x 


b = E 


xm 


where 


a 2 = 


70c^ - 12c^ + 2c^ 


16 x 


m 


-84c^ + 360^ “ 4c^ 
16x 


a 4 = 


C 1 “ 


m 


30c^ - 14 c 2 + 2c^ 


16x 


m 


- a X x m 


m 


C 2 — b 3^ 
2x b 

, m 

and c 3 = _____ 


(27) 


Now, in the region of body normal coordinates, Equations (6a, b) and 
(7) become 


u = — <J> + (1 + F)cos0 

v = gd> - (1 + F)sin0 
n 


(28a) 

(28b) 
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(1 a 2 } h ( h Ve a 2 H 8 + (1 " a 2 )g(8lJ, n ) n 


/2uv k , sin0 N f , J. r/1 U N K J. cose,^ 

+ ( T~ H + T _) H +5 + 1(1 —> H + “r ]8<() r 1 

a a 


+ sin 0 - (1 - ^-)cos0] F 

a or ** 


- [H- cose + (1 -^)sin0]gF = 0 (29) 

Z Z X] 

a a 

Likewise in the region of sheared cylindrical coordinates, Equations 
(8a, b) and (9) are transformed into the following: 


u=l + F+ f<f>^- r^g <f>^ 


(30a) 


v = g* 


(30b) 


(i -2j)f(f* ? ) s - 2f *t r i» + ¥*?„ 


+ K 1 - 4> + <r b )2(1 - V + ¥ r iie(g* n ) n 

a a a 


+ I? - r Sd - 4>ls*n + (1 -4 ,F 5 

a a 


' l r i (1 - + ^ 8F n = 0 < 31 > 

a a 

If r' = 0 and r“ = 0 in the region of body normal coordinates and 0 
b b 

= 0 and k = 0 in the region of sheared cylindrical coordinates. 
Equations (28a, b), (29), (30a, b) and (31) can be combined into a 


single set of equations. 


12 


and 


u = H + (1 + F)cos0 " r b 8<|, n 


v » - (1 + F)sin0 


' 3 . S. 




(32a) 

(32b) 


♦ K 1 -4> + (r i >2(1 - 4 + ¥ r ; ig<8 Vn 

a a a 


,2uv < , sin0. f . rn _ U \ < + cos 8 

+ <— h + ~T~ ) h *s [( ~ H + 

a a 


- r”(l - ^)]g^ + sine + (1 - ^)cos6] | 


- [HI cose + (1 - ^sine + r^(l - \)]gF n = 0 


(33) 


2.6 Botated Finite Difference Scheme 


Rotated difference (Ref. 12) is needed to keep diagonal 
dominance of the tridiagonal implicit scheme and the correct zone of 
dependence and thus the numerical stability. In the case of body 
normal coordinates, Equation (28a, b) and (29), the streamwise and 
normal derivatives <f> gs and ^ are given by 


+ss-VlT <1 V? + ¥*5n +v2 *nn 1 

q 


^NN " 2 [ H ( H Vs H V + U W 

q 


(34a) 


(34b) 
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where S and N are the streamwise and normal directions to a 
streamline. 

In the sheared cylindrical coordinates. Equations (30a, b) and 
(31), <J>g S and ^ are given as 


*SS ■ ^2 I u2 +55 + 2u(T * r i u) *5n + <v ‘ r i u)2 *nn 1 

+NN - h [v2 *55 ' 2v(u + r i v> hn + (u + r J. v) \r, 1 


(35a) 

(35b) 


Now Equation (33) is written in the form: 


n _ , . /2uv < sin9 l 

^ 2^SS + ^NN + ^ 2 H + r ^ H ^ 


+ k 1 -4 l + £ r i - r i: (1 

a a 


+ [~ sine + (1 - ■H-Ocose] F 

a a ** 

2 2 

- cose + (1 - — j)sin6 + r^(l - iL.)]gF = 0 (36) 

a a a 11 

At supersonic points, upwind differences are used for the three 

second derivatives contributing to <|>gg, and central differences are 

used for those contributing to (j^ and all first derivatives. At 

subsonic points, the usual procedure is used with central 

differences for all derivatives directly in Equation (33). Thus at 

subsonic points the truncation error is formally of the second 

order, while at supersonic points it is of the first order. 
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Equation (36) is seen to be quite similar to that used in 
RAXBOD (Ref. 12), except for the rotation function derivatives. 
Therefore, Keller and South's transonic disturbance potential code, 
RAXBOD, is modified to solve the present problem. 

2.7 Boundary Conditions 

At infinity, the perturbation potential is required to vanish; 
that is, 

f + 0 as n + » (37) 

In sheared cylindrical coordinates, the perturbations at downstream 
infinity (5 ®) must likewise vanish. This can be accomplished via 

transformation x = t( £) by mapping % - 00 to a finite value of x, or 
one can simply use a sufficiently large £ and apply <p = 0 there, or 
extrapolate <j> when M(r)^ n ^ > 1. The latter course was taken in the 
present study. That is, for M(r)^ n ^ < 1, the downstream boundary is 
located about three-fourths to one body length beyond the sting/body 
junction or other most downstream obstacle. For M ( r )^ n f > the 
only requirement is that the boundary must be downstream of the last 
subsonic region. Numerical results are otherwise insensitive to the 
precise location of the boundary. 

On the surface, n = 0, the flow tangency condition depends on 
the coordinate system as follows: 

Body Coordinates : 

v = 0 (38a) 

or 

4> - (1 + F)sin0 (38b) 
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Sheared Cylindrical Coordinates: 


v - ur' = 0 (39a) 

b 

or 

r h 

K = 5 7 (1 + F + <(, ) (39b) 

n 1 + (r£) Z ^ 

In the sheared cylindrical coordinates, the body surface 
boundary condition is satisfied by introducing dummy points inside 
the body. Details can be found in Reference 1. For completeness, 
the formulation is described in the following. Note that the dummy 
points may be located above or below the symmetry axis. For dummy 
points above the axis, as shown in Figure 2(a), the values of the 
potential function at these dummy points are computed through 
Equation (39b). 


or 


i _ ^i,jmax-l , jmax+1 

*y Ttiy 


. , - <f> /2Ay 

i,jmax+l i,jmax-l y 


(40) 


Note that = <j)^y ^ through Equation (22). It is possible that 
dummy points may be below the axis, as shown in Figures (2b) and 
(2c). Due to symmetry, the potential at a point below the axis 
should be the same as that for a point (i.e., the image point) at an 
equal distance above the axis. In this case, let y^ be the 
computational coordinate at the image dummy point where the 
potential is to be calculated. A Taylor series expansion for <{> at 
this point (which is the same as <|>^ j mav+1 ) yields (Ref. 1) 


^i, jmax+1 ^i,jmax + ^l^y + 2 ^yy + 


(41a) 
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Similarly, 


<p. . + Ay<J> + — ? — <b + 

Y i,jmax-1 Y i,jmax Y y 2 Y yy 


(41b) 


Eliminating <pyy from these equations and solving for <p^ it: 

is obtained that 


$ 


i , jmax+1 


(Ay)‘ 


<p. . , + (1 

Y i,jmax-1 


(Ay) 


tt) <>. + y, (1 

2 r i,jmax 


Ay^y 

(42) 


To calculate y^, the computational coordinate corresponding to the 
location of the axis, y a , is first obtained. Then y^ = Ay + 2y a « 
Note that y a is negative. y a can be found from the stretching 
function (Equation 22) by expansion in a series for a small y to 
give 

£ = y + ay 2 + q(0t 2 + 1} y 3 + ±- ■ 1 >< q + 2 > y 4 + . . . (43) 

Equation (43) can be inverted to give 

„ _ n _ , h\ 2 a(3a - 1) ,ru3 a(16a 2 - 12a + 2) ,r\4 

y ” A a( A ; + 2 6 ( F + * * * 

(44) 

Putting q = -r^ into Equation (44) gives the value of y a » 


2.8 Grid Halving 

A considerable saving in computing time can be achieved by 
first obtaining the solution on a coarse grid and then halving the 
mesh size in both directions for further calculation. This process 
can be continued to any desired mesh refinement within the computer 
time and storage limitations. The following third-order 
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interpolation formulas are used to interpolate results in a coarser 
grid to those in a finer grid: 


1) For points next to symmetry axis, 

♦i,j = * 5625<(> i-l,j + * 5<f, i+l,j " *° 625<(, i+3,j 
if the symmetry axis is at i - 1; 

♦i,J ' - 5625 +i+l,j + ‘ •° 625 *l-3,j 

if the symmetry axis is at i + 1. 

2) For points not next to symmetry axis, 

*i,j = ,3125 *i-i,j + * 9375 *1+1, j " * 3125< ( ) i + 3, j 

+ .0625<j, 1+5>j 

if the symmetry axis is at i - 1; 

*i,j = * 3125 *i+l,j + * 9375<f> i-l,j ” ,3125 *i-3 t j 
+ . 0625 4>. _ . 

9 J 

if the symmetry axis is at i + 1. 

Similar formulas are also used in the j direction. 

2.9 Optimization formulation: 


CONMIN (Ref. 13) is used to couple the present program 

designing an axisymmetric body. 

The objective function OBJ is formulated as 

OBJ * -0.1/(0.001 + C^) 

where C d is the pressure drag, 
w 


(45a) 


(45b) 


(46a) 


(46b) 


for 


(47) 
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The maximum thickness is assumed to be constrained. It is 


formulated as 

G( 1) = io ( Iggg_ - i) (48a) 

upper 

G(2) = 10(1 - ^ ma - X - ) (48b) 

lower 

where is the maximum radius. G is the constraint function. 

Since equality constraints are not practical numerically for 
nonlinear problems, an upper limit r upper and a lower limit r lower 
are used instead. The constant, 10, is used to increase the 
relative importance of constraint gradients in finding the optimal 
direction during optimization. 

The trailing-edge thickness can also be constrained. The 


constraint functions are defined as 

G(3) - t ^ - 1 (49a) 

upper 

G(4) = 1 - r<: - e - (49b) 

L lower 

if the constrained thickness, t, is not zero. Otherwise, they are 

G(3) - 100(r te - t upper > (50a) 

G(4) - 100(t lower - r te ) (50b) 


Since transonic computation is very CPU-intensive, the 
following representation of body shapes is used to reduce the number 
of design variables. For an ellipsoid-type body, the slope of the 
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body shape is expressed in a series as follows: 

dr A n+1 _ 0 A n+2 „ 9 . ? . . , 

j — cot j - — j — tan y + Z A^sin n0 (51a, 

X — x 

x = x^ + —-- 2 — - (1 - cos0) (51b^ 

0 0 

where cot ^ an< d tan *y ta ^ e care of the leading edge and trailing 
edge slopes, respectively; x^ is the starting x coordinate; x^ is 
the ending x coordinate; 0 is the corresponding angle in the 
transformed plane; and N is the number of coefficients in the sine 
series. The body shape can be integrated to give 
r dr , 

r • / ar dr 

“ 7 ( 8 + »in9) - -jr^ ( 9 - sin 9) + -ji ( 9 - ° 1 ^ 26 ) 


+ j, ^ .sin(n - 1)0 sin(n -I- 1) 9 . ■■ 

_ n n - 1 n + 1 ^ ! 

n=2 


By defining the following quantities 


o ~l/2x , v 

0 = cos (- r - - 1) 


x^ * y (1 + cos > M even, 1 < y < M - 1 

Weber (Ref. 14) showed that the leading edge radius r^ e is given by 

r „ M-l sin* r 

/ 2 -i£ = -2 Z (-1) P 4 - T-- ? • - (53) 

£ , ' 1 + cos* £ v ' 

y=l r ]i 

r ' M-l sin* r 

A = / 2 -iS - -2 Z (-l) y — .. H , . -H (54) 

N+l £ . ' 1 + cos* £ v ' 

y~i p 


A M4-9 = / 2 -^ = -2 Z (-l) w , Ji 

N+2 £ - 1 - cos* £ 

y=l r y 
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and 


4r 


te 


A 1 tt£ ^+1 + \+2 
Let f ( 0) = — - ( 9 + sin0) + ^±2- ( 6 - sine) 


“1 sin2 0 

— (9 —2 — ) 


^ A r sin(n - 1)0 sin(n + 1)0, 
« n l n - 1 n + 1 J 


Multiply (57) by 


sin(m -1)0 _ sin(m +1)0 
m — 1 m + 1 

and integrate with respect to 0 to obtain 

7 r 

I i 

o 


r sin(m - 1)0 sin(m + 1)0, JO 

i f(0)[ — — — ST1 — ]de 


A - A 0 A - A 
tt r m m-2 m m+2 -> 

9 2 + 2 ^ 


(m - 1)‘ 


(m + 1)‘ 


(56) 


(57) 


Hr. r 1 J. 1 

■y (A [ y + 

Z “ (m - 1 V (m + 1) 


2 ] 


H o 

m-2 


m+2 


(m - 1) (m + 1)' 


(58) 


Therefore, 


A [■ 
m l 


(m - 1) (m + l) 4 


m-2 


m+2 


(m - 1) (m + l) 4 


_ M sin(m - 1)0, sin(m + 1)0. 

+ |^ 1 f (Vt — m+l- - H > m > 2 < 59 > 

k=l 

Note that r te may become negative during the optimization 
process. In this case, the coefficient A N+ 2 is slightly reduced 
repeatedly until r te is nonnegative. 
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To find r , dr/d0 is set to zero: 


AAA 
^ = | {-j£L (1 + cose) - (1 - cos 0) + Y- (1 - cos2 0) 


N 


+ I A [cos(n - 1) 0 - cos(n + 1)0]} 


n 


= | (Vl (1 + COBS) - V2 


2 (1 - cose) + 2 ~ (1 “ cos2 0) 


N 

+ 2 E A sin n0sin0} 
2 ° 


= 0 


(60) 


Let 


T/ A N+1 ... .. \+2 

1( 0) = —7y — (1 + cose) - 


2 (1 - cos0) + 2 ~ (1 ~ cos2 0) 


Then 


I’(0) 


N 

+ 2 E A sin n0sin0 
2 n 


^N+l ^N+2 

— sin0 - — ^ — sin0 + A^sin2 0 


(61) 


N 

+ 2 E A [n cos n0sin0 + sin n0cos0] 
~ n 


(62) 


0 at r__„ can be iteratively solved by Newton's method as 

UlaA 


or 


1(0) + I'(0.)(0 1+1 


I(0 i } 

0 i+l 9 i ” I' ( 0. ) 


9 i ) - 0 , 


(63) 


For a body with nonblunt trailing edge, the shape function 
given by Equation (52) is found inappropriate because the 
coefficient A N+ ^ is too dominating. Any change in A N+ ^ will affect 
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not only the nose shape but also the trailing-edge thickness quite 
significantly. Therefore, the shape function is redefined to be 

t fl N 

r = j {A N+1 sin0cos j + £ A n cos(n - 1)6} (64) 

n=l 

It can be shown that Aj^ is still related to the leading-edge 
radius through Equation (54). A n , n < N, are determined as Fourier 
coefficients in the usual manner. 
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3. NUMERICAL RESULTS AND DISCUSSIONS 


Examination of the governing equation (Equation 9) for the 
present nonuniform flow problem indicates that the equation is 
similar to that for the uniform flow except the "nonhomogeneous” F- 
terms. Therefore, it is appropriate and convenient to modify the 
uniform-flow code of Reference 1 to solve the present problem. 

Before numerical results are presented, first some 
considerations of numerical stability and convergence of the revised 
code will be given. Relaxation and supersonic damping factors, as 
discussed in Reference 12, are needed to ensure stability and 
convergence. Therefore, they will be considered next. 

3.1 Numerical Stability 

Stability is indicated by M . Since the governing equation 

involves a i-term, A*^ occurs usually on the symmetry axis or on 

the body surface. With the addition of and F^ terms, the 

location of Ad> moves to where the maximum values of these terms 
max 

occur. The latter are somewhere ahead of the nose and away from the 
axis. When shock is strong, the solution and may be 

oscillatory. 

3.2 Convergence 

Residual of the governing equation is indicative of how well 
the current values of <{> satisfy the governing equation and thus is 
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used as the convergence criterion of the present method. In 
subsonic or low transonic free stream, the value of the residual can 
be reduced to an arbitrarily small value. However, because of the 
first-order accuracy inherent at the supersonic points, this seldom 
can be done in high transonic or supersonic freestream. The 
location of maximum residual usually occurs at either the trailing 
edge or the nose stagnation point. 

3.3 Relaxation and Supersonic Damping Factors 

A relaxation factor is used to control the stability and 

convergence at subsonic points, while a supersonic damping factor is 

to increase the stability at supersonic points. When the sum of 

residuals of the last ten iterations increases, the original code 

will increase the value of the supersonic damping factor by 0.1 or 

decrease the relaxation factor by ten percent. It turns out in most 

cases that the maximum residual occurs at either the nose or the 

tail where the flow is usually subsonic. Therefore, the supersonic 

damping factor will not change during the iteration. Since Ad> is 

max 

also an important indicator for stability and convergence and its 

location is usually not at the body surface, another indicator is 

set up to indicate whether the point with Ad> is subsonic or 

max 

supersonic. Therefore, for each ten iterations, if either the 
A^ma v~P°i n t or P°^ nt °f maximum residual is supersonic, the 
supersonic damping factor is increased by 0.1. Likewise, when 
either point is subsonic, the relaxation factor is decreased by ten 
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percent. The maximum supersonic damping factor is set to 3.0 and 

the minimum relaxation factor is set to 0.3. If for a continuous 

one hundred iterations the sum of maximum residual decreases, the 

supersonic damping factor is decreased by 0.1 when at either 

locations of maximum residual or Ad> the flow is supersonic. The 

max 

relaxation factor is increased by ten percent when at either 

locations of maximum residual or A<b the flow is subsonic. The 

max 

minimum of supersonic damping is set to zero, while the maximum 
relaxation factor is an input quantity. 

3.4 Numerical Results in Analysis 

Experimental data for a body in axisymmetric nonuniform 
transonic flow are not available for comparison. Therefore, in the 
following only theoretical results will be presented to show the 
general trend. In uniform flow, some results with data comparison 
can be found in Reference 12. 

The main motivation of this research is to find the 
nonuniformity effects on a propfan nacelle. The experimental Mach 
number profile of a propfan is plotted in Figure 3 with a scaled 
ellipsoid. Calculated pressure distributions shown in Figure 4 
indicate that the pressure distribution in a nonuniform flow is more 
negative than that in a uniform flow with a Mach number equal to 
that either of the external flow or in the slipstream. Similar 
results have been obtained for a Joukowsky airfoil in two- 
dimensional incompressible flow in Reference 16. For an ellipsoid 
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with sting as shown in Figure 5, a similar trend in pressure 
distribution as presented in Figure 6 is observed. Physically, it 
is possible that this is due to the constraint effect of the outer 
subsonic freestream which reflects the disturbance back to the 
central region. The effect of sting on the ellipsoid is similar to 
having a thick wake and is to decrease the pressure as shown in 
Figure 7. Notice that in all cases shown above, no local supersonic 
regions are present for the configuration used with a fineness ratio 
of 10. 

In Reference 8, Mach numbers of 0.98 and 1.1 were used in 
determining axisymmetric bodies with minimum pressure drag in 
uniform flow. Therefore, nonuniform transonic freestreams from Mach 
0.98 to 1.1, 1.2, 1.3, and 1.4 are chosen in the present parametric 
investigation. The following Mach profiles will be used (see 
Figures 8a, b): 

M^ n ^(r) = 1.4 - 0.42 tanh(^) (62) 

where r is the radial distance and d controls the extent of the 
nonuniformity region. As shown in Figure 9, it is seen that for the 
same maximum Mach difference, the pressure distribution appears to 
be about the same, irrespective of difference in profile shapes. 

Note that Cp in Figure 9 and all that follow is based on the dynamic 
pressure in the external uniform flow. This result is unexpected 
because in References 15 and 16 the extent of nonuniformity was 
shown to affect the pressure distribution of an airfoil in two- 
dimensional flow. To investigate this problem further, step-type 
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nonuniform profiles shown in Figure 10 are employed. Again, the 
same results are obtained as shown in Figure 11. 

On the other hand, for the configuration with sting in the same 
Mach profiles (Figure 12), some differences (Figure 13) do show 
up. However, for the step-type nonuniform profiles (Figure 14), all 
pressure distributions, again, are the same (Figure 15). In Figure 
16, the sting effect is seen to reduce the shock peak. It also 
shows that the nonuniform Mach profile shape is not an important 
parameters in the present problem. One possible reason for this is 
that there is no vortex flow in the present axisymmetric cases, i.e. 
with zero lift, while in the airfoil problem (Refs. 15 and 16), lift 
is significant. 

Since the nonuniform Mach profile shape is not an important 
parameter in the present study, nonuniform extent of d = 0.1 will be 
used in the following to investigate the effect of nonuniformity 
magnitude. As shown in Figure 17, in supersonic nonuniform 
freestreams the magnitude of nonuniformity increases the pressure 
coefficient negatively and nonlinearly for the ellipsoid 
configuration. Similar trend can be seen for the ellipsoid/sting 
configuration. In a transonic nonuniform freestream, increasing the 
Mach number in the nonuniform region tends to make C p more negative 
and move the shock rearward as shown in Figures 19 and 20. In 
Figure 21 it can be observed that the pressure is more negative in a 
subsonic outer stream and that the sting reduces the shock peaks. 
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In Figure 22, the drag coefficient is plotted for all cases 
investigated. It is seen that it is slightly negative for near- 
sonic nonuniform cases. This is because of the neglect of viscous 
drag and base drag. Transonic freestream is found to induce lower 
drag for near-sonic cases but higher drag for stronger 
nonuniformities. This is because at near sonic conditions transonic 
freestream wave drag is minimal. However with higher Mach numbers 
in the nonuniform region, the wave drag approaches that of a uniform 
supersonic freestream. 

3.5 Numerical Results in Design Optimization 

Chan's (Ref. 8) numerical results were obtained at Mach 0.98 
and 1.1 in uniform flow. However, the starting shape to be used in 
the present investigation has been verified experimentally (Ref. 17) 
and numerically (Ref. 12) not to induce a shock wave until M = 

0.986. Therefore, it is decided to use Mach 0.995 and 1.1 as 
typical Mach numbers for uniform subsonic and supersonic cases, 
respectively, and a nonuniform freestream varying from a Mach number 
of 0.995 to 1.1 for the nonuniform flow case. To reduce the number 
of design variables, representation of the body shape in a Fourier 
series as discussed in Section 2.9 is used. The number of design 
variables (i.e. the Fourier coefficients) can be reduced to six 
without sacrificing the accuracy. Since in transonic computation it 
takes too many iterations to converge when a variable is perturbed, 
the step size can not be too large. However, if a small step size 
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is used, the gradient of the objective function would be small so 
that the objective function will change little. Therefore, a user f s 
judgment is needed in the design process. The following results are 
obtained after many cycles of optimization. In each cycle, only one 
to three iterations in CONMIN are used. The results may not 
represent the final optimum. 

Ellipsoid 

For the uniform supersonic case of Mach 1.1, the original and 
the final pressure distributions are compared in Figure 23. The 
original shape produces a gradual expansion until a tail shock is 
encountered. The designed shape results in a wavy pressure 
distribution ending with sudden expansion and shock at the tail. A 
drag reduction of 14 percent is achieved. The shapes are compared 
in Figure 24.. The designed shape shows a two-percent reduction in 
maximum thickness, with thickness reduced in the front, and the 
contour straightened in the rear. 

For uniform subsonic freestream of Mach 0.995, a drag reduction 
of 5% and a minor pressure change, caused by reduction in shock 
strength, are seen in Figure 25. The designed shape shows slight 
thinning in the front and slight thickening in the rear, as shown in 
Figure 26. 

For the case with nonuniform transonic freestream, it is 
observed in Figure 27 that the design pressure distribution 
eliminates the double shocks associated with the original shape. 

This results in a 13% reduction of drag and 1.18% increase of 
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maximum thickness. In Figure 28, the designed shape is shown to be 
slightly thinner in the front, thicker in the middle, and flatter in 
the rear. 

In all cases, the location of maximum thickness stays the same 
in the subsonic case and shifts slightly forward in supersonic and 
transonic cases. Note that the original shape is shockless until M 
= 0.986. All cases considered here involve higher Mach numbers. It 
appears that by thinning the front and thickening the rear (to 
reduce the surface slope), the pressure drag can be reduced at the 
higher Mach numbers. 

Body with NACA 0012-Type Contour 

To design a body with a rounded leading edge and a trailing 
edge which is not blunt, an initial shape given by the NACA-0012 
airfoil contour is chosen. Again, six design variables (A R in 
Equation 64) are used. The maximum thickness is constrained to be 
between 13% and 11.5% and the trailing-edge thickness between 0% and 
1%. As indicated in Reference 11, reducing the residual of the 
governing equation to a small value may not be needed for a 
reasonably accurate solution. Therefore, in the following design 
process, the convergence criterion is based on the maximum equation 
residual obtained in the analysis of the input shape. The final 
designed shape is then subject to further analysis through 1200 
iterations for final plotting. 

In a uniform flow with M = 0.98, the results are presented in 
Figures 29 and 30. As can be seen in Figure 29, the shape of the 
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NACA 0012 contour produces higher negative pressure behind the nose 
and a stronger shock which is more forward. By decreasing the nose 
radius and increasing the thickness in the aft portion (Figure 30), 
the designed shape produces less expansion to reduce the negative 
pressure level behind the nose and a weaker shock which is more 
aft. The achieved pressure drag reduction is about 46% with Ac^ = 
0.0187. 

The same initial shape is again used in a nonuniform transonic 
flow. The Mach number in the external stream is 0.98. However, 
over an extent of nonuniform flow region equal to one-half of the 
body length, M is set to 0.995 around the body. Similar results in 
pressure distributions and change in body shape are obtained, as 
shown in Figures 31 and 32. That is, reducing the nose radius and 
increasing the thickness in the aft portion tend to reduce the wave 
drag. However, the pressure drag reduction is less than that in the 
uniform flow case, being 29% with Ac d « 0.0137. 
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4. CONCLUSIONS 


An inviscid nonuniform transonic axisymmetric body code capable 
of performing analysis and design was developed. Numerical stabil- 
ity and convergence behaviors were discussed, and so were the super- 
sonic damping and relaxation factors. Numerical results showed that 
nonuniformity caused pressure coefficient to be more negative. 

Sting attached to the body was to reduce the pressure peak near the 
juncture. If a shock was present, the strength was reduced and its 
location moved forward. The extent and shape of the nonuniformity 
region appeared to have little effect on pressure distribution. 
Increase in nonuniformity magnitude would make Cp more negative and 
the shock location more rearward. 

The CONMIN optimizer was coupled with the present analysis code 
to design axisymmetric bodies in uniform and nonuniform flow. For 
an ellipsoid, the trend indicated that by thinning the front portion 
and flattening the rear of a body, the pressure drag could be re- 
duced at high transonic and low supersonic speeds. The drag reduc- 
tion in a uniform flow of M = 1.1 and 0.995 was 14 percent and 5 
percent, respectively. In a nonuniform flow of M * 0.995 to 1.1, 
the pressure drag reduction achieved was 13 percent. For a body 
with a rounded leading edge and nonblunt trailing edge, the nose 
radius should be reduced and the thickness in the aft portion 
slightly increased to decrease the pressure drag. Using the NACA 
0012 contour as the initial shape, it was shown that a drag reduc- 
tion of 46 percent and 29 percent was achieved, respectively, in a 
uniform flow of M = 0.98 and a nonuniform flow of M = 0.98 to 0.995. 
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Figure 2a. Ordinary Dummy Point 
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Figure 2b. Image Dummy Point inside the Body 
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Figure 2c. Image Dummy Point outside the Body 



Mach Number 

Figure 3. Propfon Mach Number Profile and the 
the Relative Location of an Ellipsoid 
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Figure 4. Propfon Slipstreom Effect on 

Pressure Distribution of on Ellipsoid 
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Figure 5. Propfan Mach Number Profile and the 
Relative Location of an Qllpsold/Stlng 
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Figure 6. Propfan Slipstream Effect on Pressure 
Distribution of an Ellipsoid/Sting 
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Figure 7. Comparison of Pressure Distributions 

between Ellipsoids with and without Sting 
In a Propfon Stream Given in Figure 5 
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Figure 8a. Mach Number Profiles for Different Spread Parameters 
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Figure 8b. Mach Number Profiles for Different Spread Parameters 
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Figure 9. Spread Effect of Transonic Nonuniform Freestreams 

on Pressure Distributions of an Ellipsoid 
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Figure 10. Mach Number Profiles for Different Nonuniformity Thickness 


47 



- 1.5 



x/c 

Figure 11. Depth Effect of Transonic Nonuniform Freestreoms 

on Pressure Distributions of an Ellipsoid 
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Figure 12. Mach Number Profiles for Different Spread Parameters 
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Figure 13. Spread Effect of Transonic Nonuniform Freestreams 

on Pressure Distribution of an Ellipsoid/Sting 
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Figure 14. Mach Number Profiles for Different Nonuniformity Thickness 
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Figure 15. Depth Effect of Transonic Nonuniform Freestreams 

on Pressure Distribution of an Ellipsoid/Sting 


52 



Pressure Coefficient Cp 



«/c 

Figure 16. Comparison of Pressure Distributions 
between Ellipsoids with/without a Sting 
ond Different Types of Non uniformities 
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Figure 17. Effect of Supersonic Nonuniform Freestreams 

on Pressure Distribution of an Ellipsoid 
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Figure 18. Effect of Supersonic Nonuniform Freestreoms 
on Pressure Distribution of an Ellipsoid/Sting 
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Figure 19. Effect of Transonic Nonuniform Freestreoms 

on Pressure Distribution of on Ellipsoid 
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Figure 20. Effect of Transonic Nonuniform Freestream on 
Pressure Distribution of an Ellipsoid/Sting 
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Figure 21. Comparison of Pressure Distributions 
of Ellipsoids wfth/wfthout Sting In 
Supersonic and Transonic Free Stream 
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Figure 22. Effect of Nonuniformity Magnitude on Drag 

Coefficient. Symbols at Higher Mach Numbers 
Represent the Condition of Uniform Flow. 
Symbols at Higher Mach Numbers Represent 
Nonuniform Flow With Maximum Mach Numbers at 
the Indicated Values. 
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Figure 23. Comparison of Original and Design Shapes and 
Pressure Distributions in Uniform Row of M-1.1 
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Figure 25. Comparison of Original and Design Shapes and 
Pressure Distributions In Uniform Flow of M**.995 
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Comporiaon of Original and Design Shapes and Pressure 
Distributions in a Transonic Nonuniform Freestream 
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Figure 29. Comparison of Original and Design Shapes and 
Pressure Distributions In Uniform Flow of M=0.98 
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Figure 31. Comparison of Original and Design Shapes and Pressure 
Distributions In a Transonic Nonunrform Freestream 
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